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Abstract 

We have demonstrated that rate dependent restitution and action potential duration-refractory 
period hysteresis can be reproduced in a one-dimensional two-variable Chernyak-Starobin-Cohen 
reaction-diffusion medium with variable excitation threshold. We show that restitution and hys- 
teresis depend on the relationship between pacing period and steady state excitation threshold and 
also on the rate of excitation threshold adaptation after an abrupt change in pacing period. It was 
also observed that the onset of action potential duration alternans is determined by the minimal 
stable wavefront speed, which could be approximated by the analytical critical speed of a stable 
solitary pulse. This approximation was suitably accurate regardless of the adaptation constant of 
excitation threshold, its dependence on pacing interval, or magnitude of the slopes of restitution 
curves. 
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I. INTRODUCTION 



Repetitive pacing of biological reaction diffusion media by over-threshold stimuli gives 
rise to periodically propagating excitation waves. In cardiac tissue the duration of excitation 
referred to as action potential duration, Tap, as well as the speeds of excitation wavefront and 
waveback, depend on previous stimulation periods and refractory (diastolic) intervals, To/ 
[l, 2, Is]. The analysis of such dependences known as restitution curves has been established 



as an effective method for evaluation of normal functioning of the heart 

Pioneering experimental studies 1^, [ll, 12] established two major pacing sequences such 
as steady state (dynamic) and S1-S2 pacing protocols, which led to two standard Tap resti- 
tution dependences. It has been shown that the protocol dependent rates of Tap adaptation 
were different for stepwise S1-S2 perturbations of cycle length acceleration and deceleration 



121 . Il3| . This phenomenon was later experimentally introduced as action potential dura- 



tion cycle length hysteresis |l4j. Ilq and has been recently associated with cardiac ischemia, 
coronary flow reduction 16|, [l7|, and cardiac memory ISl. Il9]. 

A stepwise change in pacing rate following a long series of conditioning SI stimuli may 
result in prolonged adaptation of action potential duration to its new steady-state value. 
The process of such adaptation can extend well beyond the Tap response to the first S2 test 



stimulus 
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12| and give rise to the additional constant BCL 



restitution component attributed to cardiac memory 



Dasic cycle length) transient 
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221 ]. Relationship between 



this phenomenon and stability of pulse propagation as well as with rate dependent Tap-Toi 



lysteresis has 
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3een recently investigated in experimental 
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231 ] and theoretical 



28| studies. 



It has been found that in the presence of restitution transients cardiac dynamics is more 
complex than predicted by Nolasco and Dahlen restitution criterion 29[ . Specifically, it was 



demonstrated that dynamic restitution curve slopes > 1 may not automatical 



of excitation wave stability and subsequent appearance of alternans 20 



24 



ly indicate loss 
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27 



30]. 



Computational experiments with Fenton-Karma and Mitchell-Schaeffer ionic models readily 



identified that such effects can be quantified by fitting T, 



curves using time dependent gate variables 
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and Tpti for specific restitution 



33]. 



In this paper, we analyze excitation wave propagation in a one- dimensional cable based 
on the approach, which follows from direct experimental observations of the dependence of 
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cardiac muscle resting potential on frequency of external pacing [3J, l35|, [SGj. We implement 



a two-variable exactly solvable Chernyak-Starobin-Cohen (CSC) react ion- diffusion model 



37 



and modify it accordingly to incorporate pacing rate driven adjustments of resting 
potential as a rate dependent excitation threshold. We use an exponential-like evolution of 
excitation threshold, V^, that takes place over the course of multiple heart beats following 



stepwise changes in pacing rate j35|. We demonstrate that adaptation of K significantly 
affects the stability of pulse propagation and gives rise to rate dependent Tap-Tdi hysteresis. 
In particular, we find that under given medium parameters, regardless of the slopes of 
restitution curves, the appearance of alternans is determined by proximity of the wavefront 



speed to the minimal speed of a stationary solitary pulse determined analytically in |37l |. 



II. METHODS 

Basic equations that describe a class of exactly solvable models for excitable media have 



been defined in 



37 



38| . Here we introduce a modification of this analytical model by 
adjusting the excitation threshold, V^, in response to changes in frequency of external pacing. 
We will consider the model in dimensionless form: 

^ = + (1) 



for u < V 
ioi u> V 



^ = e{Cu + Vr-v) (2) 

dVr _ -Vr + B{t) 
Itt ~ T 

u{x,t) and v{x,t) are a membrane potential and slow recovery current, respectively. A, e, (, 
and T are the model parameters, where < e. The scaling of the system is described in 
the Appendix. 

The pacing function, P{x,t), is defined as a product of two functions, X{x) and T(t). 
Each is composed of the Heaviside step function, 9(x), as follows: X{x) = A[Q{x — 6i) — 
6(x — ^2)] and T(t) = 6(tA;) — 6(tfc+T5), where A and 52 — ^1 are the amplitude and width of 
the pulse, respectively, Tg is the pulse duration, and = tAr(m.-i) + [k — N{m — l)]Tm are the 
instants of time when stimuli are delivered. We use this to construct a pacing protocol where 
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the pacing period, Tm, is stepwise constant (Tm > Tg). N represents the number of stimuh 
at each pacing interval plateau, the index m denotes each pacing plateau, m = 1, . . . , M, 
and k is an integer in the range k = 1, . . . , NM. The number of stimuli N is the same 
for all plateaus. Overall during the course of the protocol, progressively decreases to a 
minimum and then increases to its starting value. 

The right-hand term B in Eq. [3] responds to the stepwise evolution of pacing period, Tm- 
Stepwise changes in B result in smooth exponential transition of the excitation threshold 
from one steady-state plateau to another. A steady-state value of excitation threshold at 
each m^^ plateau, VJ^ = B^, was chosen to be linearly dependent on the corresponding 



pacing interval, |34j. |35|. 



5± = -/3±T^ + a± (4) 

Here and are positive parameters that determine the amplitude of change of B^ 
between two consecutive pacing plateaus for the increasing and decreasing rate, respectively. 

III. NUMERICAL SIMULATIONS 

The system of Eqs. [T]|3] was solved numerically on a short cable of 250 grid points with 
spatial and temporal grid intervals of Ax = 0.13 and At = 7.2 x 10~^, respectively. The 
length of the cable was approximately equal to the width of the pulse to reflect the relative 
dimensions of the heart and a propagating cardiac pulse at moderate heart rates. Periodic 
wavetrains were produced by stimulating the cable with a square wave at the left end using 
the function P{x,t) = X{x)T{t), defined above, where A = 10, 6i = 2Ax, 62 = 15Ax, 
and Ts = 10'^ At. The model parameters A, e, and ( were equal to 0.4, 0.1, and 1.2, re- 
spectively, for all simulations. Numerical solutions were computed using a second-order 



explicit-difference scheme jo]. A typical solution is depicted in Fig. [H showing the propaga- 
tion of a single pulse at three instants of time. 

In order to quantify the dynamics of the system (Eqs. [I]l3]), we computed the action 
potential duration. Tap, the diastolic interval, Tq/, and the wavefront velocity, c. The 
action potential duration was defined as the interval of time when m > f at a specified 
node, Xq. Accordingly, the refractory period, Tp/, was defined as the interval of time when 
u < V. The speed of the wavefront, c, was calculated based on the time it took for a point of 
constant phase on the wavefront {u = 0.5) to travel a span of 10 grid points centered around 
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FIG. 1: Spatio-temporal dynamics of u,v variables. Three snapshots of u and v are shown for a 
excitation threshold Vr = 0.19. The first snapshot occurs at the instant of the pacing stimulus and 
shows the formation of the wavefront. The next two snapshots are taken after periods of 7.2 (20% 
of the pacing period) and 14.4 and show formation and propagation of the pulse. 

xq. In order to analyze a developed pulse, we measured intervals and speeds at xq = 20 
(except Sec. IIIIDI) where the speed of the wavefront had reached a constant steady-state 
value. 



A. Restitution for constant excitation threshold 



The system of equations [T]l3] was initially studied with Eq. [3] replaced by its asymptotic 
form ^ = 0, which is equivalent to r = oo. The cable was stimulated periodically for 
forty consecutive pacing periods over the range = 70 to 25 with decrements of 1.5. At 
the end of each plateau, Tap and To/ had reached steady-state values that were used to 
compose the steady-state restitution curve. Two values of Tm, at 30 and 27, were used to 



obtain conventional S1-S2 restitution curves 



10 



ll| . The pulse, T^p , resulting from the 



test stimulus, S2, following the conditioning sequence and the last diastolic interval, T^j, 
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FIG. 2: Steady-state restitution curves are shown with sohd hues, and S1-S2 restitution curves 
computed for SI BCLs of 27 and 30 are shown with dashed hues. Critical Tap and Tdi values at 
which pulse durations begin to oscillate are at the ends of the curves where T^"* = 19.1, T^'p = 5.7 
and r57* = 20.8, Tj/ = 5.5 for Vr = 0.19 and K- = 0.215, respectively. Open markers in the 
upper left insert indicate corresponding critical speeds plotted with the analytical dispersion curve 
for a solitary pulse [37|]. The mid and lower inserts demonstrate oscillation of pulse duration at 



and below the critical speed. 



from the conditioning SI plateau composed the S1-S2 restitution curve. 

Steady state and S1-S2 restitution curves computed for a pair of constant excitation 
thresholds and a pair of SI basic cycle lengths are shown in Fig. [2l We observed that the 
maximal difference between steady state and S1-S2 restitution curves comprised less than 
15% of the corresponding steady state value of Tap, which implied that transient responses 
to any premature stimulus were on average limited to just a single non-stationary pulse. 
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The ends of both steady state restitution curves shown in Fig. [2] indicate the critical points 
below which no stable propagation and no 1 : 1 responses were observed. We found that 
when the pacing interval reached the value of = 26.3 at = 0.215, the resulting action 
potentials oscillated in duration as shown in the mid level insert panel. Further reduction of 
the pacing interval from = 26.3 to = 25.9 induced an even more complex 3 : 2 T^p 
response pattern depicted in the lower level insert panel. 

The upper level insert shows the dispersion curve computed analytically for a steady state 
solitary pulse [37|. The critical speed and duration of a stable solitary pulse for the model 
parameters, A, (, and e, described above is Ccrit = 0.48 and T^p = 5.4. As shown in the 
insert, the difference between the numerical and analytical critical speeds is relatively small 
and amounts only to 10% of the wavefront speed determined from our numerical model at 
Vr = 0.215. This suggests a criterion for determining the critical speed below which pulse 
durations start to oscillate. On the contrary, the conventional Nolasco-Dahlen critical slope 
stability criterion is unsuitable as the maximum slope of the dynamic restitution curve at 
this Vr is 30% greater than one 



B. Restitution and hysteresis for rate-dependent Vr 

Unlike the previous section, the evolution of Vr according to Eq. [3] resulted in a set of 
prolonged transients initiated by abrupt changes in stimulation rate. These Tap transients, 
which constituted the constant BCL restitution (negative slope. Fig. [3]A), had a duration 
of 5-50 stimulation periods depending on the adaptation constant, r (Fig. [3p). The phase 
of constant BCL adaptation followed the immediate S1-S2 responses, which were on the 
contrary positive and aligned with the steady state restitution curves for each steady state 
excitation threshold (gray lines. Fig. [3]A.). 

Different steady state values of excitation threshold during progressively increasing and 
decreasing stimulation rates gave rise to Tap-Tdj hysteresis. The cable was stimulated 
periodically for fifty consecutive periods at a series of SI conditioning plateaus. The series 
consisted of five plateaus with decreasing followed by the same number of plateaus with 
increasing Tm- The excitation threshold evolved according to Eq. [HI and its steady state 
value at each plateau was related to the pacing period by Eq. |H 

Using this protocol we demonstrated that higher values of Tap-Tdj hysteresis corre- 
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FIG. 3: Restitution relations and Tap-Tdi hysteresis for a stepwise stimulation protocol with 
rate dependent excitation threshold. (A) A series of five stepwise changes of pacing rate are 
shown with individual Tap and Tdj as black circles. Open circles mark steady-state Tap and Tdj 
intervals. Portions of constant excitation threshold, = B^, and steady state restitution curves 
are shown in gray. (B) Tap-Tdi hysteresis produced by increasing differences between and 
B^. /3+ = 6 X 10~^ and = 0.37 for the top curve (solid line, squares), /3~ = 4 x 10~^ and 
a" = 0.31 (dotted line, solid circles), P~ = 2 x 10~^ and a~ = 0.25 (dashed line, triangles), and 
P" = 0.1 X 10~^ and a" = 0.19 for the lowest curve (dash-dotted line, circles). (C) Tap-T^i 
hysteresis for different values of the adaptation constant, r (r = 216, dotted line; r = 32, solid 
line), when /3+ = jd" and a"*" = a~. (D) Adaptation of action potential duration after an abrupt 
decrease in stimulation interval for two different rate constants {f3~^ = 6 x 10~^ and a"*" = 0.37). 
The measured transient intervals in Panel C are highlighted with corresponding markers. 
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sponded to greater differences between and during the decremental and incremental 
stages of the pacing protocol (Fig. [3]B). We also found that Tap-Tdi hysteresis was depen- 
dent on the adaptation constant, r, if Tap and To/ were measured before their steady-state 
values were reached. For instance, for a time lag of seven stimulation intervals (Fig. EP), 
the magnitude of Tap-T^u hysteresis increased as the adaptation constant increased from 
r = 32 to r = 216 (Fig.[3p). 



C. Transitions to Tap alternans, and comparison with stability of a solitary pulse 

When the pulse duration and speed are less than certain critical values, steady-state 



solitary pulse and wavetrain propagation in an infinite cable do not exist 37], |38|. Such 



critical values for a solitary pulse occur at the end of the solitary pulse's dispersion curve. 



as computed analytically in 37|]. Similar critical values for our short cable were described 
earlier in Sec. IIII Al for a medium with constant excitation threshold. In this section, we 
analyze perturbations of Tap and wavefront speed when steady state pulse speeds are below 
a certain value close to the critical speed of a solitary pulse. 

We perturbed both pulse duration and wavefront speed near their analytical critical 
values, T^p and Ccrit, as indicated by arrows labeled "A,B" and "C,D" (Fig.Hj). If the speeds 
of the wavefronts elicited after perturbation exceeded Ccrit less than 22%, we observed Tap 
alternans (inserts A and B) that developed between two values indicated by open markers at 
the end of the arrow "A,B" . Alternans occurred sooner when the adaptation rate was higher 
due to a smaller adaptation constant, r (insert A). The minimal non-alternating wavefront 
speeds computed at xq = 20 for different were the same regardless of changes of 5^ or 
the four fold increase of the steady-state restitution slope (Fig. [5]). On the contrary, when 
the perturbation resulted in a wavefront speed that exceeded Ccrit by more than 22% (the 
end of the arrow "C,D"), similar bifurcations did not occur. Instead, we observed typical 
exponential adaptation from one steady-state to another (inserts C and D). 



D. Influence of propagation on the spatial distribution of hysteresis and alternans 

We demonstrate that Tap adaptation to a stepwise change in pacing interval is different 
for different points of observation along the cable. When the distance between the point 
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FIG. 4: Solitary pulse dispersion curve and initiation of T^p alternans at xq = 20. Inserts (A) 
and (B) show Tap response to a perturbation in stimulation rate, for initial = 46.8, final 
Tm = 40.3, initial = 0.31, and final = 0.32. The adaptation constant is small, r = 32, 
in A and large, r = 216, in B. Inserts (C) and (D) depict stable transitions of Tap for the same 
adaptation constants (r = 216 in C; r = 32 in D) when the starting and ending pulse speeds 
are greater than in (A) and (B) (initial = 51.3, final = 45.0, initial B^ = 0.31, and final 
= 0.32). Gray line in main panel depicts dispersion curve of a steady-state solitary pulse [371]. 
The markers at the left end of the dispersion curve show the starting speed and Tap (solid circles) 
and the ending speeds and Tap (open circles) for both inserts. 

of observation and stimulation site increases, the wavefront and waveback speeds decrease. 
At the point Xq = 6 adjacent to the stimulation site, for an initial = 46.8 and final 
Tm = 40.2, the wavefront speed, c = 0.69, is substantially higher than the critical speed of 
a solitary pulse, Ccrit = 0.48. Under these conditions. Tap gradually adapts to a new steady 
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FIG. 5: Steady state restitution for three different sets of rate dependence parameters, (3~^ and . 
T = 216 for all curves, xq = 20. The ends of the curves indicate the values of T^p below which 
the duration of pulses oscillates. Note that the speed of the pulses at the ends of the curves is the 
same regardless of the magnitudes of restitution slopes, which vary with parameters (3^ and a^. 

state value as shown in Fig. [6]A.. On the contrary, at xq = 20, for the same change in pacing 
interval and excitation threshold parameter, the wavefront speed is substantially lower, 
c = 0.59, which results in a series of oscillating Tap- The lower branch of alternating Tap 
(Fig. [6]B) corresponds to slowly propagating wavefronts whose speeds are virtually equal 
(7% difference) to the analytical value of Ccrit- Further increase of the distance between 
the observation point and stimulation site results in the increase of the amplitude of Tap 
alternans. The closest point at which alternans can be observed is located at the midpoint 
of the cable (Fig.ET). 

We also observed that the magnitude of Tap-Tdj hysteresis is larger when measured 
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FIG. 6: The spatial distribution of 1 : 1 and 2 : 2 responses along the cable. Oscillations were 
induced by perturbing Tap, with initial = 46.8, final Tm = 40.2, initial Bm = 0.31, and final 
Bm = 0.32. Panels A and B show the response to perturbation at = 6 and xq = 20, respectively. 
Panel C shows the spatial distribution of the last two Tap of the train of 50 pulses illustrated in 
Panels A and B. 

further from the stimulation site. Figure [7] shows the twofold increase of the magnitude of 
hysteresis and Tap between two observation points at Xq = 4 and Xq = 20. 



IV. CONCLUSIONS 



We have demonstrated that rate dependent restitution and Tap-Tdi interval hysteresis 
can be reproduced in a one-dimensional two- variable CSC reaction-diffusion medium where 
excitation threshold adjusts to changes in pacing rate. We show that the rate dependence of 
restitution and hysteresis are influenced by two major factors. The first one is the adaptation 
constant, r, of excitation threshold evolution after an abrupt change in pacing interval. The 
second is the dependence, B^, of the steady-state excitation threshold on the pacing period. 
We show that steady-state and S1-S2 restitution curves coincide if the excitation threshold 
is constant, which corroborates with earlier findings for other reaction diffusion models 

B, y, 



28l |. On the contrary, the steady-state and S1-S2 restitution curves diverge if the 



magnitude of excitation threshold varies with changes in pacing rate, leading to prolonged 
Tap transients following an abrupt change in pacing interval. 

We found that larger values of Tap-Tdi interval hysteresis were associated with greater 
differences between 5+ and B^. Even if there was no difference between these values, 
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FIG. 7: Steady-state Tap-Tdj hysteresis at xq = 4 (dashed hne, open circles) and xq = 20 (sohd 
hne, squares). For both curves, /3+ = 6 x 10"'^, a+ = 0.37, /3~ = 2 x 10~^, and a' = 0.25. 

hysteresis resulted from increasing adaptation constants, r, if Tap and Ta were measured 
before their steady-state values were reached. 

Our numerical simulations showed that stimulating the cable with short pacing intervals 
and high excitation thresholds elicited slower pulses that led to Tap alternans. It was 
observed that the minimal stable wavefront speed could be approximated by the analytical 



critical speed of a stable solitary pulse 37|, ISSj. This approximation was suitably accurate 
regardless of values of and magnitudes of the slopes of restitution curves. We also found 
that the onset of alternans occurring after an abrupt change in pacing rate was more delayed 
for larger values of r. 



APPENDIX: 



The scale of u is the maximum steady-state action potential amplitude Uq, the scale of 
V is given by cTfUo, and the time scale is Cm/ erf, where af corresponds to the maximum 
sodium conductance and Cm is the membrane capacitance. The characteristic length scale is 
given by y^UJoy, where D is the diffusion coefficient. The small parameter e <C 1, is equal 
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to Cm/ (TsCf) and C = CTs/af where as corresponds to the maximum potassium conductance. 
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